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Often gaining insight into the functioning of biomolecular systems requires to follow their dynamics along 
a microscopic reaction coordinate (RC) on a macroscopic time scale, which is beyond the reach of current all 
atom molecular dynamics (MD) simulations. A practical approach to this inherently multiscale problem is to 
model the system as a fictitious overdamped Brownian particle that diffuses along the RC in the presence of 
an effective potential of mean force (PMF) due to the rest of the system. By employing the recently proposed 
FR method [I. Kosztin et al., J. of Chem. Phys. 124, 064106 (2006)], which requires only a small number of 
fast nonequilibrium MD simulations of the system in both forward and time reversed directions along the RC, 
we reconstruct the PMF: (1) of deca-alanine as a function of its end-to-end distance, and (2) that guides the 
motion of potassium ions through the gramicidin A channel. In both cases the computed PMFs are found to be 
in good agreement with previous results obtained by different methods. Our approach appears to be about one 
order of magnitude faster than the other PMF calculation methods and, in addition, it also provides the position 
dependent diffusion coefficient along the RC. Thus, the obtained PMF and diffusion coefficient can be used in 
a suitable stochastic model to estimate important characteristics of the studied systems, e.g., the mean folding 
time of the stretched deca-alanine and the mean diffusion time of the potassium ion through gramicidin A. 



PACS numbers: 87.15.A-, 87.10.Tf, 87.10.Mn, 05.70.Ln, 05.40.Jc 



I. INTRODUCTION 

The study of the structure-function relationship of 
biomolecular systems often requires to follow their dynamics 
with almost atomic spatial resolution on a macroscopic time 
scale, which is beyond the reach of current all atom molecular 
dynamics (MD) simulations. A typical example is molecu- 
lar and ion transport through channel proteins [1]. Indeed, 
in order to determine the forces that guide the diffusion of 
molecules across the channel one needs to know with atomic 
precision the structure of the channel protein-lipid-solvent en- 
vironment. However, the duration of the permeation process 
across the channel occurs on a time scale (e.g., /xs to ms) that 
may exceed by several orders of magnitude the time scale of 
several tens of nanoseconds currently attainable by all atom 
molecular dynamics (MD) simulations[2]. Whenever the dy- 
namic properties of interest of such a system can be described 
in terms of a small number of reaction coordinates (RCs) then 
a practical approach to this inherently multiscale problem is 
to model the system as fictitious overdamped Brownian parti- 
cles that diffuse along the RCs in the presence of an effective 
potential of mean force (PMF) that describes their interaction 
with the rest of the system. 

Recently we have proposed an efficient method for calculat- 
ing simultaneously both the PMF, U{R), and the correspond- 
ing diffusion coefficient, D, along a RC, R, by employing a 
small number of fast nonequilibrium MD simulations in both 
forward (F) and time reversed (R) directions |3J. The ef- 
ficiency of this method, refeiTed to as the FR method, was 
demonstrated by calculating the PMF and the diffusion coef- 
ficient of single-file water molecules in single-walled carbon 
nanotubes ^. The obtained results were found to be in very 
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good agreement with the results from other PMF calculation 
methods, e.g., umbrella sampling |4, 5, 6]. 

To further test its viability, in this paper we apply the FR 
method to investigate the energetics of two well-studied ex- 
emplary systems, i.e., (i) the helix-to-coil transition of deca- 
alanine in vacuum, and (ii) the transport of K^ ions in the 
gramicidin A (gA) channel protein, inserted in a fully sol- 
vated POPE lipid bilayer In each case we seek to calculate 
the PMF as a function of a proper RC, i.e., the end-to-end 
distance {R) of deca-alanine and the position (z-coordinate) of 
the potassium ion along the axis of the gA channel. The com- 
puted PMFs are found to be in good agreement with previous 
results obtained by using either the Jarzynski equality ||7|,l8l] or 
the umbrella sampling method fllsla]- However, compared 
to these PMF calculation methods our approach is about one 
order of magnitude faster and, in addition, also provides the 
position dependent diffusion coefficient along the RC. Thus, 
by employing the computed PMF and diffusion coefficient in 
a suitable stochastic model we could estimate important char- 
acteristics of the studied systems, e.g., the mean folding time 
of the stretched deca-alanine and the mean first passage time 
of K^ through the gA channel. 

The remaining of the paper is organized as follows. To 
make the presentation self-contained, in Sec.|ll]a brief descrip- 
tion of the FR method is provided, along with the theory used 
to analyse our results. The study of deca-alanine is described 
in Sec. |III1 while that of K^ transport in the gA channel in 
Sec.|lV] Finally, Sec.|V]is reserved for conclusions. 



II. THEORY 

By definition, for a classical mechanical system described 
by the Hamiltonian Hq{T), the PMF (Landau free energy), 
U{R), along a properly chosen RC (R) is determined from the 
equilibrium distribution function of the system by integrating 



out all degrees of freedom except R, i.e., ij] 

. „-PHoir) 
e-P"^^^ ^ po{R) = /^r— 5[R-R{r)] . (!) 



Here poiR) is the equilibrium distribution function of the RC, 
Zq is the partition function, j3 = l/ksT is the usual thermal 
factor, and 8{R) is the Dirac-delta function whose filtering 
property guarantees that the integrand in Eq.([T]) is nonzero 
only when R{r) — R. In this paper we use the convention that 
R [or R{t)] is the target value, while R = R{r) is the actual 
value of the RC. Also, it is convenient to use kgT as energy 
unit. Thus, in Eq. ([TJ one needs to set j3 = 1 . 

Unfortunately, by using equilibrium MD simulations the 
direct application of Eq. ^ is practical only for calculating 
U{R) about its local minimum. An efficient way to properly 
sample R is provided by steered molecular dynamics (SMD) 
[9] in which the system is guided, according to a predefined 
protocol, along the RC by using, e.g., a harmonic guiding po- 
tential 



Vr{R{T)) = ''-[R{T)-R]\ 



(2) 



where k is the elastic constant of the harmonic guiding poten- 
tial. With this extra potential energy, the Hamiltonian of the 
new biased system becomes Hr = Hq +Vr{R). As a result, 
atom "/' in the selection that define the reaction coordinate 
will experience an additional force 



dVR 



-k[R{T)-R] 



dm 



(3) 



By choosing a sufficiently large value for the elastic constant 
k, i.e., the so-called stiff-spring approximation QIISIj the dis- 
tance between the target and actual value of the RC at a given 
time can be kept below a desired value. 

In constant velocity SMD simulations |9], starting from an 
equilibrium state characterized by R{0), the target value of the 
RC (or control parameter) R{t) is varied in time according to 
R{t) ~ R{Q) + vt, < t < T, where v is the constant pulling 
speed. For each such forward (F) path there is a time reversed 
(R) one in which the system starts from an equilibrium state 
corresponding to R{t) and reaches ^(0) according to the pro- 
tocol RR{t) = /Jf (t - f) = R{t) - vf, < f < T. The external 
work done during a SMD simulation is given by 



Wf= dR [dVR{R)/dR] = k / dR{R - R 



R(t) 
Ro 



(4) 



The F and R work distributions are not independent but related 
through the Crooks Fluctuation Theorem 11 1111 



PfiW) 



„w., 



Pr{-W) 
where the F dissipative work is given by 

W^f^Wf-AU , 



(5) 



(6) 



with AU = U[R{t)] - U[R{0)]. In principle, the PMF can be 
determined from the so-called Jarzynski equality (JE) I,1Z1 



(exp(-Wdf )) = 1 , 



(7) 



that follows directly from Eq. ^ 1 1 1]. Within the stiff-spring 
approximation the sought PMF is given by the second cumu- 
lant approximation ||3, IM llOll 



AUpiR) 



-log{exp[-WFiR)])^{WF) 



-aj/2. 



(8) 



o} = {W^)-{WFf 



where CTp is the variance (2nd cumulant) of the F work. Also, 
within the stiff-spring approximation the work distribution 
function Pf{W) is Gaussian and, therefore, the cumulant ap- 
proximation (|8]l is exact ||7|l. However, in practice Eq. (O is 
valid only close to equilibrium because SMD pulling paths 
can sample only a narrow region about the peak of the Gaus- 
sian Pf {W), while the validity of JE is crucially dependent on 
very rare trajectories with negative dissipative work (Wd < 0). 
Thus, in general, having only a few SMD trajectories one can 
determine fairly accurately the mean work (Wf) but not the 
variance Gp, which in most cases is seriously underestimated. 
In the FR method this shortcoming is eliminated by com- 
bining both F and R pulling trajectories and employing 
Eq. (|5]l, which is more general than the JE (|7]i. Within the stiff- 
spring approximation, Eq. (|5]l implies that the F and R work 
distribution functions are identical but displaced Gaussians, 
and the PMF and the mean dissipative work Wj = Wjf = WdR 
can be determined from the following simple equations |3;] 



AU^i{WF)-{WR))/2 



{W,i) = {{Wf) + {Wr))/2. 



and 



--2{Wd). 



(9a) 



(9b) 



(9c) 



Equations (|9]l are the key formulas of our FR method for cal- 
culating PMFs from fast F and R SMD pullings. Clearly, the 
superiority of the FR method, for calculating the PMF (and the 
mean dissipative work), compared to the one based on the JE 
equation is due to the fact that Eqs.(|9|l contain only the mean F 
and R work (whose values can be estimated rather accurately 
even from a few SMD trajectories) and not the corresponding 
variance. In fact the latter (see Eq. ( |9c] l) is also determined by 
the mean F and R work. 

Although, strictly speaking, the FR method can only deter- 
mine the PMF difference between initially equilibrated states 
connected by F and R SMD trajectories, in practice we find 
that in many cases Eqs (|9a]i-(|9b]i give good results even be- 
tween the division points 7?,, / = 1, . . . ,A^— 1, of the interested 
interval [Rq = R{0),Rn = R{t)]. The reason for this is that 
for a stiff harmonic guiding potential the equilibrium distri- 
bution of the RC is a narrow Gaussian that can be sampled 
through very short MD simulations. Thus, even if the system 
is far from equilibrium due to fast pulling by a sufficiently 
stiff spring, the instantaneous value of the RC will always be 



sufficiently close to its equilibrium value. However, even in 
such cases the pulling speed should not exceed values that 
would cause excessive perturbation to the rest of the degrees 
offreedomof the system. Thus, the number of division points, 
A^, does not need to be large, implying a fairly small com- 
putational overhead for the equilibration of the system at Ri, 
i = l,...,N-l. 

An alternative approach for calculating the PMF difference 
between two equilibrium states connected by np forward and 
riK reverse SMD paths is based on the maximum likelihood es- 
timator (MLE) method applied to Crooks' fluctuation theorem 
© fij, i.e.. 



£t 



1 



J i I «/r/«Rexp(Wf, — At/) 



(10) 



St 



+ n/;/«/7 exp(- Wft- - At/) 



= 0. 



We use Eq. ( fTOl i to test the accuracy of the PMF results ob- 
tained with our PR method. 

Finally, since it is reasonable to assume that W^i is propor- 
tional to the pulling speed v, one can readily determine the 
position dependent friction coefficient y{R) from the slope of 
the mean dissipative work y{R) = {d{W(j{R))/dR) /v. Then, 
the corresponding diffusion coefficient is given by the Einstein 
relation (in ksT energy units) ^ 



D{R) = 7(/?)-' = v{d{W,iiR))/dR) 



-1 



(11) 



Once both U{R) and D{R) are determined, the dynamics of 
the reaction coordinate on a macroscopic time scale can be 
described by the Langevin equation corresponding to an over- 
damped Brownian particle (ij] 



Y{R)R = -dU{R)/dR + £,{t), 



(12a) 



or equivalently, by the corresponding Fokker-Planck equation 
for the probability distribution function p{R,t) of the reaction 
coordinate 



d,p{R,t) = -dRJ{R,t) 



dR[D{R)dRp{R,t)]+dR[U'{R)p{R,t)] , 



(12b) 



where £,{t) is the Langevin force (modeled as a Gaussian 
white noise) and j{R,t) is the probability current density. 

For example, Eq. ( I12bl ) can be used to calculate the mean 
folding time of deca-alanine (see Sec. IIII Bb from the com- 
pletely stretched (coil) conformation 7?^ = 33 A to the folded 
(helical) conformation Rq = 14.5 A as the corresponding 
mean first passage time (MFPT) 11511 . i.e.. 



T= rdRe''^''^/D{R) TdR'e"^'"'^ 



Re 



R, 



(13) 



III. STRETCHING DECA-ALANINE 

Deca-alanine is a small oligopeptide composed of ten ala- 
nine residues (Fig.[T]). The equilibrium conformation of deca- 
alanine, in the absence of solvent and coupled to an artificial 
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FIG. 1 : (color online), (a) Cartoon representation of deca-alanine. 
The reaction coordinate R is defined as the distance between the 
first (CAj) and last (CAio) Ca atoms, i.e., the end-to-end distance 
of the peptide. The spring, with elastic constant k, connecting 
CAi and CAio corresponds to an elastic guiding potential V(R;t) = 
{k/2) [R — Roit)] that can be used to cycle deca-alanine between the 
(b) folded and (c) unfolded (completely stretched) conformations. 
In (b) and (c) the backbone (sidechain) atoms are shown in cartoon 
(CPK) representation. In the folded (b) configuration the hydrogen 
bonds that stabilize the a-helix are also shown. (Snapshots rendered 
with the program VMD lISll ). 



heat bath at room temperature, is an a— helix. The system 
can be stretched to an extended (coil) conformation by ap- 
plying an external force that puHs its ends apart. Once the 
stretched system is released it will refold spontaneously into 
its native a— helical conformation. Thus, this can be regarded 
as a simple protein unfolding and refolding problem that can 
be comfortably studied via SMD simulations due to the rel- 
atively small (104 atoms) system size. It is natural to de- 
fine the reaction coordinate as the distance R between the first 
(CAi) and the last (CAio) Ca atoms. To calculate the PMF, 
U{R), that describes the energetics of the folding/unfolding 
process we have use SMD simulations to generate a small 
number (in general 10) F and R pulling trajectories and ap- 
ply the PMF calculation methods described in Sec.|lll i.e., the 
FR method [Eqs. ©], the JE method [Eq. ^] and the MLE 
method [Eq. (fTotl. The SMD harmonic guiding potential ^ 
coiTesponded to an ideal spring of tunable undeformed length 
R{t) inserted between CAi and CAio (see Fig.[T^). Note that 
this choice of the guiding potential is more natural than the 
one customarily used in the literature in which the atom at- 
tached to one of the two ends of the spring is fixed [|8l [l6l[l7ll . 



A. Computer modeling and SMD simulations 

The computer model of deca-alanine was built by em- 
ploying the molecular modeling software VMD |18]. All 
simulations were performed with NAMD 2.5 119] and the 
CHARMM27 force field for proteins lH^, SH]. A cutoff of 
I2A (switching function starting at lOA) for van der Waals 
interactions were used. An integration time step of 2 fs 
was employed by using the SHAKE constraint on all hy- 
drogen atoms [2^. The temperature was kept constant (at 
300 K) by coupling the system to a Langevin heat bath. The 
system was subjected to several equilibrium MD and non- 
equilibrium SMD simulations. We divided the reaction coor- 



dinate R G [13,33] A into ten equidistant intervals (windows) 
delimited by the points 7?,- = (13 + 2/) A, / = 0, . . . , 10. Next, 
a pool of equilibrium states were generated for each Rj from 
4 ns long equilibrium MD trajectories. These states were used 
as starting configurations for the SMD F and R pulls on each 
of the ten intervals. The spring constant in these equilibrium 
MD simulations was ^ — 50 kcal/mol/A^. The equilibrium 
length of the folded deca-alanine was determined from two 
free MD simulations starting from a compressed (R — 13 A) 
and the completely stretched (R — 33 A) configurations of 
deca-alanine. Both simulations led to the same equilibrium 
length Rec, = 14.5 A. 

In order to calculate U{R) a total of six sets of F and R 
SMD simulations were carried out. In each of the first three 
sets of SMD runs we used ten simulation windows, but three 
different pulling speeds: vi = 1 A/ps, V2 = 10^' A/ps, and 
vo = lO^'* A/ps. The sets corresponding to vi,2 consisted of 
10 F and 10 R SMD trajectories. For the quasi-equilibrium 
pulling speed v'o only one F and R runs were performed. In the 
last three sets of SMD simulations we used a single simulation 
window, covering the entire range of the RC, and used the 
same three pulling speeds as in the previous SMD runs. For 
all six sets of SMD simulations, the stiff-spring constant was 
k = 500 kcal/mol/A2. 

To construct the forward and reverse work distribution 
functions on the segment R e [17,21] A, we performed 2000 
F and the same number of R SMD simulations. In order to 
generate a sufficient number of starting equilibrium config- 
urations it was necessary to extend the equilibration runs at 
both R= n A and 7? = 21 A to 5 ns. In all these simulations 
we used a pulling speed of v = 1 A/ps and a spring constant 
ofA: = 500kcal/mol/A2. 

Finally, to estimate the mean refolding time of the com- 
pletely stretched deca-alanine we performed 100 free MD 
simulations starting from an equilibrium configuration corre- 
sponding to /? = 33 A. As soon as deca-alanine reached its 
folded, equilibrium length Rgq = 14.5 A the simulation was 
stopped and the refolding time recorded. 



B. Results and Discussion 

The PMFs calculated using the FR method corresponding 
to the six different pulling protocols described in Sec. 1111 Al 
are shown in Fig. |2] As expected, for the very small pulling 
speed V'o the system is in quasi-equilibrium throughout the 
SMD runs leading to the same (true) PMF regardless of the 
number of simulation windows considered. However, while 
the dissipative work is negligible for both F and R processes, 
repetition of these simulations resulted in different PMFs for 
R>24- A, and it will be discussed below (see also Fig.|5]l. Not 
surprisingly, in case of the very fast pulling speed vj , the PMF 
for the single simulation window is rather poor along R ex- 
cept at the end-points of the window. Indeed, the FR method 
allows to calculate the PMF difference between two equilib- 
rium states connected by fast F and R SMD processes that 
follow the same protocol. However, it is remarkable that us- 
ing ten simulation windows, even at this large pulling speed 






50 
40 
'30 
20 
10 



T 1 \ \ \ \ 1 \ \ \ 

v=1 0"" A/ps [1 segs; 1 pull] 

7=10"" A/ps [1 seg; 1 pull] 

v=10"' A/ps 110 pulls/seg/dir] 

v=10"' A/ps[1 seg; 10pulls/dir] 

. — . v=1 A/ps [10 pulls/seg/dir] 
v=1 A/ps [1 seg; 10 pulls/dir] 

A v=10"' A/ps[MLE] 

o v=1 A/ps [MLE] 




J I I I I I I L. 



R[A] 

FIG. 2: Potential of mean force (PMF) of deca-alanine as a function 
of the reaction coordinate R. The different curves were obtained with 
the FR method by employing different simulation and PMF calcula- 
tion protocols described in the text. 



the resulting PMF is rather close to the real one. For the still 
fast pulling speed V2 the situation is similar. While the single 
simulation window case lead to a rather poor PMF (though 
somewhat better than in the v'l case), the ten simulation win- 
dows result is almost indistinguishable from the true PMF. 
For comparison, the PMFs calculated at /?,, / = 1 , . . . , 10 using 
the MLE method for both vj and V2 are also shown in Fig. |2] 
Based on these results one may conclude that the FR method 
gives very good PMF even for fast pulling speeds and using 
only a few F and R trajectories, provided that a sufficient num- 
ber of simulation windows are used. 

A comparison between U{R) obtained from the FR method 
and the cumulant approximation of the JE method (applied 
separately for the F and for the R SMD trajectories) are shown 
in Fig. [3] In general, the FR method yields better PMF in all 
cases, and especially when one employs (i) one simulation 
window (Fig. [3^ and c), and (ii) a very large pulling speed vj 
(Fig.[3]a and b). For the ten simulation windows with pulling 
speed V2 (Fig. [3]d) the FR and JE methods are comparable 
though even in this case the JE F (R) method systematically 
over (under) estimates the PMF. Note, however, that an av- 
erage of the JE PMFs for the F and R trajectories leads to a 
result very close to the FR one. 

An important prediction of the FR method is that, provided 
that the stiff-spring approximation holds, the F and R work 
distributions are identical Gaussians centered about the mean 
F and R work, and therefore shifted by 2 At/. To test this pre- 
diction we have determined the work distribution histogram 
corresponding to 2000 F and a same number of R SMD trajec- 
tories corresponding to the RC segment R £ [17,21] A. The 
results are shown in Fig. H) Although the histograms seem 
to be Gaussian (dashed lines) they are not identical as pre- 
dicted by the FR method. In a previous study 1 17] the clear 
deviation from Gaussian of the external work distribution in 
case of deca-alanine was pointed out and it was attributed to 
the non-Markovian nature of the underlying dynamics of the 
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FIG. 3: Comparison between the PMFs U{R) obtained using the 
PR method (thick solid line) and the cumulant approximation of the 
JE corresponding to the ten forward (dashed line) and reverse (dot- 
dashed line) SMD trajectories, respectively. The thin solid line cor- 
responds to the exact PMF. The upper (lower) panels correspond to 
a uniform pulling speed of 1 A/ns (0. 1 A/ns). The PMFs in the right 
panels were determined by dividing the 20 A pulling distance into 
ten equidistant segments (the system being equilibrated in each of 
the end points of the individual segments), while the PMFs in the left 
panels were determined by considering the entire pulling distance as 
a single segment. 



system. However, in our case both work distributions look 
Gaussian and the relatively small but clearly noticeable dif- 
ference between them may be due either to the failure of the 
stiff-spring approximation or to incomplete sampling. After 
all, the end-to-end distance is a poor and insufficient reaction 
coordinate for describing the folding and unfolding processes 
of a polypeptide. 

This last point becomes rather clear when the system is sub- 
jected to repeated folding (R) and unfolding (F) processes at 
the quasi-equilibrium speed vq. At this speed the system is 
at almost equilibrium throughout the SMD pulls and one ex- 
pect that the PMF is given by the external work, i.e., the dis- 
sipated energy (which is a stochastic quantity) is negligible. 
While for R <2A A one gets systematically the same PMF, 
for R G [24, 33] A one obtains different PMFs depending on 
the direction of pulling, as one can see in Fig. |5j). A care- 
ful inspection of these trajectories reveal that the folding and 
unfolding processes occur through different pathways in the 
above mentioned range of the RC. Thus, it appears that R is 
not sufficient to specify the metastable intermediate states of 
the system, and a more complete description requires the in- 
troduction of extra order parameters, e.g., the distribution of 
the hydrogen bonds (H-bonds) in the peptide. Indeed, the dy- 
namics of the formation and rupture of the H-bonds during 
folding and unfolding, respectively, may be rather different. 
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FIG. 4: Histogram of the distribution functions (thin solid lines) of 
the forward (Wf) and reverse (W«) works along the segment R e 
[17, 21] A. Although the histograms seem to be Gaussian (dashed 
lines) they are not identical as predicted by the FR method (see text 
for details). 
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FIG. 5: (a) Variation of the number of hydrogen bonds in deca- 
alanine during the quasi-equilibrium (v = 10^ A/ns) F and R 
pullings. (b) The PMF U(R) calculated as the external work done 
during the quasi-equilibrium F (dashed line) and R (solid line) 
pullings. The discrepancy between the two PMFs is most likely due 
to the difference on how the H-bonds are formed and destroyed dur- 
ing the forced folding and unfolding processes, respectively, as indi- 
cated in the inset snapshots of the peptide. Dark (lite) color corre- 
sponds to the R (F) process. 



As shown in the inset snapshots in Fig. [Sj), the formation of 
the six H-bonds during the R process is much more homoge- 
neous than their rupture during the corresponding F process. 
This observation is reinforced by the time dependence of the 
average number of H-bonds in deca-alanine shown in Fig.|5^. 
Thus there are at least two distinctive pathways in the helix- 
to-coil transition of deca-alanine, both being explored during 
quasi-static pullings. During fast pulling, however, one of the 
pathways is preferred compared to the other 

Finally, as an application of the determined PMF and the 
diffusion coefficient, which was found to be approximately 

constant D « 0.27 A /ps, we calculated the mean folding time 
(i.e., coil-to-helix transition) by employing Eq. ( fTSl ). The the- 
oretical result of T w 140 ps compares rather well with the 
MFPT of « 100 ps obtained from the 100 free MD refolding 
simulations described in Sec. lIIIAl 



IV. K+ TRANSPORT IN GRAMICIDIN A CHANNEL 

Gramicidin A (gA) is the smallest known ion channel that 
selectively conducts cations across lipid bilayers [i23il . gA is 
a dimer of two barrel-like j3 -helices that form a ~ 26 A long 
and 4 — 5 A wide cylindrical pore through the lipid membrane 
(Fig. |6]l. Each helix consists of 15 alternating Asp and Leu 
amino acids. Due to its structural simplicity, gA is an im- 
portant testing system for ion permeation models, and it has 
been extensively studied in the literature both experimentally 
and through computer modeling. NMR studies have shown 
that each end of the channel has a cation binding site that is 
occupied as the ion concentration is increased [24] . The con- 
ductance is at maximum when the average ion population in 
the channel is one. The backbone carbonyls inside the pore 
are oriented such that the electronegative oxygen atoms face 
inward. The cation selectivity of gA is mainly due to these 
oxygens, which attract cations and repel anions 112511261 12711 . 

In spite of its structural simplicity, the energetics of the 
ion transport through gA is far from trivial. Computation- 
ally, most of the difficulty arises from the sensitivity to errors 
due to finite-size effects and from the poor description of the 
polarization effects by the existing force-fields. Besides the 
cation gA also accommodates ^ 6 single-file water molecules 
II28II (see Fig.|6]3) whose arrangement and orientation seems to 
play an important role in stabilizing the ion within the channel 

Previous PMF calculations of the potassium ion, K^, 
through gA yielded a large central barrier that resulted in a 
conductance orders of magnitude below those measured. In 
has been speculated that the measured conductance can be re- 
produced by a PMF that has a ~ 8 k^T deep energy well at 
both ends of the channel and a ^^ 5 k^T barrier in the middle 
II30I1 . Although PMF calculation methods that try to compen- 
sate for finite-size and polarization effects have improved in 
recent years, they continue to yield results that do not match 
the experimental ones. Most of these methods employ equilib- 
rium MD simulations with umbrella sampling ll3ll[32l[33ll and 
combined MD simulations with continuum electrostatics the- 
ory 1.34.1 . Recent attempt to apply the JE method (see Sec.HIli 




FIG. 6: (color online), (a) Gramicidin-A channel (/3-helix dimer col- 
ored in green) in POPE lipid bilayer (grey), solvated in water (van 
der Waals representation). The potassium ions is shown as a blue 
sphere, (b) Cross-section of the gramicidin A channel. The A'+ ion 
(blue) and the water molecules move single file inside the pore. (Cre- 
ated with the program VMD 118,1 ). 



for calculating the PMF of A'+ in gA did not yield the desired 
result [35f]. Here we apply our FR method to calculate both 
the PMF, U{z), and the position dependent diffusion coeffi- 
cient, D{z), of K^ in gA, and compare our results with the 
ones from the literature. 



A. Computer modeling and SMD simulations 

The computer model of gA was constructed from its high 
resolution NMR structure (Protein Data Bank code UNO 
1361 ). After adding the missing hy drog ens, the structure was 
energy minimized. Using the VMD IBTT plugin Membrane the 
system was inserted into a previously pre-equilibrated patch 
of POPE lipid bilayer with size 72 x 72 A^. Lipids within 
0.55 A of the protein were removed. Then, the membrane- 
protein complex was solvated in water, using the VMD plugin 
Solvate. The final system contained a total of 36,727 atoms, 
including 155 lipid molecules and 5, 700 water molecules. Af- 
ter proper energy minimization and 0.5 ns long equilibration 



of the system, a A'+ ion was added at the entrance of the chan- 
nel. To preserve change neutrality a Cl^ counterion was also 
added to the solvent. Finally, the system was again energy 
minimized for 10,000 steps and equilibrated for 0.5 ns with 
K^ placed in three different positions along the z-axis of the 
channel, namely at z G {—15,0, 15} A. The origin of the z- 
axis corresponded to the middle of gA (see Fig. |6j5). In order 
to prevent the pore from being dragged during the SMD pulls 
of the K^ ion, two types of restraints were imposed: (i) back- 
bone atoms restrained to their equilibrium positions (referred 
to as fully restrained); and (ii) backbone atoms restrained only 
along the z-axis (referred to as z- restrained). 

The F and R SMD simulations (needed to obtain the PMF 
using the FR method) were performed on three systems: (SI) 
backbone of the channel fully restrained with only one pair of 
K^ and Cl^ ions in the system; (S2) backbone of the channel 
fully restrained with 200 mM electrolyte concentration (ob- 
tained by adding 20 extra pairs of K^ and Cl^ ions to the 
solvent using the VMD plugin Autoionize); and (S3) back- 
bone of the channel z-restrained and electrolyte concentration 
200 mM. A total of 10 F and 10 R SMD pulls were performed 
along the z-axis of gA on two segments: z G [—15,0] A and 
z G [0, 15] A, corresponding to the two helical monomers. The 
pulling speed was v = 15 A/ns, while the spring constant of 
the harmonic potential that guided K^ across the pore was 
k = 20 kcal/mol/A2. 



B. Results and discussion 

A comparison of the PMFs of K^ along the axis of gA ob- 
tained for systems SI, S2 and S3 by employing the FR method 
is shown in Fig. |7] For gA with fully restrained backbones 
(i.e., systems SI and S2) the PMFs have only a weak de- 
pendence on the electrolyte concentration, and exhibit a huge 
central potential barrier of ~ 40 kgT, which is due to the ar- 
tificially imposed rigidity of the system. Once the flexibility 
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FIG. 7: Comparison of PMFs obtained for systems SI (thin solid 
line), S2 (dotted line) and S3 for two different pulling protocols: (i) 
pulling force on K^ applied along the z-direction (dashed line), and 
(ii) K^ pulled along the axis of the channel (thick solid line). 

of the gA channel in the plane of the membrane is restored by 
restraining the backbone atoms only along the z-axis (i.e., sys- 
tem S3), the central barrier of the PMF decreases to ^ 15 kgT, 



as shown in Fig. [T] (thick solid line). The transverse flexibil- 
ity of the channel leads to fluctuations in its radius that fa- 
cilitate the diffusion of K^ along the pore. This is in total 
agreement with previously published results, which empha- 
size the crucial role played by the flexibility of the gA chan- 
nel in its cation transport properties 11321 l33l 1381] . The PMF for 
system S3 was determined with the FR method by employ- 
ing two different pulling protocols. First, the pulling force 
on K^ was applied along the z-axis (dashed line in Fig. |7]) 
but there was no restrain on the cation's motion in the cross 
section of the pore (i.e., in the xy-plane). In the second set 
of pullings, beside the elastic pulling force oriented along the 
z-axis, the potassium ion was constrained to move along the 
axis of the channel (thick solid line in Fig.|7]). As one can see 
in Fig. |7] both pulling protocols yielded essentially the same 
PMF. Thus, we preferred using routinely the second pulling 
method especially because during the first one the potassium 
ion occasionally escaped between the two helices into the lipid 
bilayer. 

The PMF, U{z), was calculated separately for the two seg- 
ments (corresponding to the two helical monomers) using 
Eqs. (|9]l. The work done during the F and R SMD pullings are 
plotted in Fig. [8^ and Fig. [8}5, respectively. Due to the sym- 
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FIG. 8: Results of the FR method calculations for system S3 and 
SMD pulling speed v — 15 A/ns. (a) individual (thin lines) and 
mean (thick line) work for F pulls; (b) individual (thin lines) and 
mean (thick line) work for R pulls; (c) U(z) along the two seg- 
ments (dashed line); the symmetrized PMF is shown as solid line 
(see text); (d) mean dissipative work VK^(z) along the two segments, 
zG (-15,0) A (dotted line) and z e (0,15) A (dashed line), and their 
arithmetic mean (solid line). The slope of the linear Wci{z) yields a 
constant diffusion coefficient D = 10.34 A^/ns. 

metry of gA with respect to its center, the PMF for the two 
segments (dashed lines in Fig. [8};) form nearly mirror-images. 
Therefore, a better estimate of the PMF for the entire gA can 
be obtained by symmetrizing U{z) with respect to the center 
of the channel (i.e. z = A) (solid line in Fig.[8];). The F and 

R mean dissipative works, W^ (z), (averaged over the two 
segments) are also shown in Fig.|8}l (dotted and dashed lines. 
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respectively). The fact that Wj (z) and Wj{z) closely match 
each other is another indication that our FR method seems to 
work fine in the case of the gA channel too. Note that Wj(z), 
averaged over the F and R processes, (thick line in Fig. |8}J) 
is almost linear, which according to Eq. (fTTT i yields a constant 
diffusion coefficient D « 10.3 A^/ns. Now, the obtained D and 
U{z) can be used to solve Eqs. (I12al i and/or ( I12bl i for making 
prediction on the long time dynamics of the K^ ion in the gA 
channel. 
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FIG. 9: Comparison between the FR method (thick solid lines) and 
the cumulant approximation of the JE approach (thick dashed lines) 
for the: (a) potential of mean force and (b) mean dissipative work, 
obtained from a small number (only 10) of F and R fast SMD pulling 
trajectories. The JE method was employed in both F (dotted lines) 
and R (dash-dotted lines) directions. 



The comparison between U{z) and Wd{z) obtained from the 
FR method (thick solid lines) and the cumulant approximation 
(CA) of the JE approach (thick dashed lines), respectively, is 
shown in Fig. |9^. The bias in the cumulant approximation 
of the JE method applied either to the F (CA/r, dotted lines) 
or to the R (CA^, dash-dotted lines) processes is manifest in 
Fig.|9] While the former (CA/r) systematically underestimates 
the peaks in the PMF and the corresponding mean dissipative 
work, the latter (CAr) systematically overestimates the same 
quantities. The difference between the central barrier height 
of the CA/r and CA/; PMFs is 3.5 k^T, while at the channel en- 
trance the difference is almost twice as big (7 kgT). The neg- 
ative (positive) bias in CA/- (CA/;) is due to the fact that the 
JE approach uses explicitly the variance (i.e., the 2nd cumu- 
lant) of the corresponding non-equilibrium work distributions, 
which (unlike the mean work) cannot be accurately estimated 
from a few SMD pullings (see Sec. |II]l. However, by averag- 
ing CA/7 and CA/; the opposite biases more or less cancel out 
and the resulting mean PMF (thick dashed line in Fig. |9^) be- 
comes a close match to U (z) calculated from the FR method. 
According to Fig. |9j3, the same conclusion can be drawn for 
the mean dissipative work as well. 



Our U{z), calculated using the FR method, (thick solid line 
in Fig. [Tol l has two ^ 6 hgT deep wells positioned at the en- 
trances in the channel iz ~ ±10. SA) and two high barriers 
of ^ 15 ksT positioned close to the center of the channel 
(z K, ±3 A). Another small barrier (~ 1.4 hgT) appears to be 
located between the two high barriers, right at the geometri- 
cal center of gA. This small center barrier is well separated 
by the two main ones by a potential well of ^^ 3.5 kgT. Ac- 
cording to Fig. [TOl our PMF (thick solid line) is rather sim- 
ilar to the ones reported in recent publications by Bastug et 
al [32] (double-dot-dashed line) and by Allen et al IBlll (dot- 
dashed line). These authors used the standard umbrella sam- 
pling (US) method f^,^^ to calculate their PMFs. As shown 
in Fig. [TOl besides the small difference in the positions of the 
wells at the ends of the gA channel, there are two notable dif- 
ferences between the PMFs obtained by the FR and US meth- 
ods. First, the barrier height of the PMF computed with the 
FR method is only ^15 k^T as compared to ^ 20 kgT ob- 
tained from US. Second, the central peak in U{z) obtained 
from the FR (US) method is ^ 2 k^T below (above) the two 
main peaks. 
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FIG. 10: PMFs of K+ in the gA channel determined from the FR 
method (by employing different SMD pulling protocols as described 
in the text), and the umbrella sampling method Blllsal . 

To test the reliability of the FR method for determining 
U{z), besides the standard pulling protocol (involving 10 
SMD pulls in both F and R directions with a pulling speed 
V = 15 A/ns), we have used two additional ones, involving 
only 5 SMD pulls in both F and R directions. The two pulling 
protocols differed only in their pulling speeds, namely v = 
15 A/ns in the first (thin-solid line in Fig.fTOb and v = 30 A/ns 
(dotted line in Fig. [TOb in the second. As seen in Fig. [TOl all 
three FR method calculations yielded a consistent PMF, with 
noticeable differences only around the ends of the gA channel. 

Although the FR method leads to U{z) similar to the US 
result (albeit with a smaller main barrier height) none of these 
PMFs is suitable for reproducing the experimentally measured 
K^ conductivity of the gA channel. This would require a 
channel entrance well depth of ~ 8 k^T and a main barrier 
height of ^ 5 k^T ll30ll . The main problems in getting these 
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values are due to the limitations of the currently used MD 
methods that use empirical non-polarizable forcefields and, 
therefore, cannot account for the induced polarization in the 
lipid hydrocarbons and, most importantly, for the polarization 
of water in the course of the MD simulations OUmill . 

In order to mimic polarization effects caused by the passage 
of K^ through the channel, we reduced the partial charge of 
the ion from +e to +0.5e in system S3 (see Sec. IIV Ab . and 
carried out new SMD F and R pullings for recalculating the 
PMF through the PR method. The resulting U{z) is shown in 
Fig. [To] (dashed line). As one can see, in the new PMF the 
potential wells at the entrance of the channel moved by 2.5 A 
towards the center and their depth increased to 8.5 hgT . Fur- 
thermore, in a more dramatic change, the height of the barrier 
decreased from ^15 k^T to ^ A.2 k^T. Although the above 
approach to account for polarization effects is rather simplis- 
tic, the obtained PMF (apart from the new positions of the 
potential wells) has the previously estimated form [ 3011 that 
is capable for describing quantitatively the transport of A'+ in 
gA. 



V. CONCLUSIONS 

In this paper we have shown that the FR method |3] pro- 
vides an effective approach for calculating both the PMF, 
U{R), and the diffusion coefficient, D{R), along a properly 
chosen reaction coordinate R, in biomolecular systems by us- 
ing only a small number of fast forward and time reversed 
constant velocity SMD simulations. The obtained PMFs for 
deca-alanine are in goo d agreement with the ones reported 
in recent studies 1^ [la] . We have found that computation- 
ally the FR method is more efficient and accurate than similar 
PMF calculation methods, e.g., the one based on the Jarzynski 
equality. By employing the computed PMF and diffusion co- 
efficient in a suitable stochastic model we could estimate im- 
portant characteristics of the studied systems, e.g., the mean 
folding time of the stretched deca-alanine. 

We also applied the FR method to calculate the PMF of 
a potassium ion through the gramicidin A channel. As ex- 
pected from previous umbrella sampling calculations, the ob- 
tained PMF featured a main central barrier of height '^ 15 k^T 
and two wells at the entrance in the channel with depth ^ 
6 ksT. The PMF was reproduced rather well when using 
a smaller number of SMD pulling trajectories and/or higher 



SMD pulling speeds, confirming the reliability of the FR 
method. The channel protein flexibility, maintained in the 
SMD simulations by restraining the corresponding backbone 
atoms only along the axis of the channel, has been shown to 
play a major role in the transport of K^ in gramicidin A. In- 
deed, the height of the main potential barrier in a rigid channel 
is almost three times higher than in the flexible one. The dis- 
sipative work inside the channel was found to be linear in z, 
yielding a constant diffusion coefficient D « 10.3 A^/ns. The 
PMF calculated from the same SMD pulls using Jarzynski's 
equality with the cumulant approximation yielded inconsis- 
tent results for both forward and reverse directions. However, 
the biases in these to directions almost cancel out when av- 
eraging the forward and reverse PMFs, leading to almost the 
same result as the FR method. Furthermore, the FR method 
yielded consistently PMFs similar to the ones using the tra- 
ditional umbrella sampling method but in considerably less 
time (i.e., - 3 days per PMF on a 64 CPU, 2.8GHz Intel Xeon 
EM64T, cluster). However, the conduction of the channel can- 
not be reproduced with any of the computed PMF profiles, 
mainly because of the very large central barrier. The main 
problem in determining PMFs in ion channels through MD 
simulations is the poor treatment of polarization effects by the 
current non-polarizable forcefields. To account for the polar- 
ization of K+ inside the channel, its effective point charge was 
reduced to +0.5e. The recalculated PMF exhibited barrier and 
well sizes very close to the values needed to reproduce the ex- 
perimental data. Hopefully, with new polarizable force fields 
the FR method will provide a simple to use, efficient and reli- 
able tool for calculating PMFs for ion and molecular transport 
through channel proteins. 
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